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Abstract 

We evaluated the performance of GBLUP including dominance genetic effect (GBLUP-D) by estimating variances and 
predicting genetic merits in a computer simulation and 2 actual traits (T4 and T5) in pigs. In simulation data, GBLUP-D 
explained more than 50% of dominance genetic variance. Moreover, GBLUP-D yielded estimated total genetic effects over 
1.2% more accurate than those yielded by GBLUP. In particular, when the dominance genetic variance was large, the 
accuracy could be substantially improved by increasing the number of markers. The dominance genetic variances in T4 and 
T5 accounted for 9.6% and 6.3% of the phenotypic variances, respectively. Estimates of such small dominance genetic 
variances contributed little to the improvement of the accuracies of estimated total genetic effects. In both simulation and 
pig data, there were nearly no differences in the estimates of additive genetic effects or their variance between GBLUP-D 
and GBLUP. Therefore, we conclude GBLUP-D is a feasible approach to improve genetic performance in crossbred 
populations with large dominance genetic variation and identify mating systems with good combining ability. 
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Introduction 

Genomic selection refers to the use of genome-wide dense single 
nucleotide polymorphism (SNP) markers to predict breeding 
values and subsequently select individuals [1]. Several approaches 
of genomic prediction have been presented. One of them is the 
genomic best linear unbiased prediction (GBLUP), which uses 
genomic information in the form of a genomic relationship matrix 
that defines the additive genetic covariance between individuals 
[2,3]. The genomic relationship coefficients are estimated with 
higher accuracy than when using pedigree information because 
genomic information can capture of Mendelian sampling across 
the genome. GBLUP has become popular approach in genomic 
selection of dairy cattle [4,5] because it is simple and has low 
computational requirements [6,7]. 

Most published models only include additive genetic effects [8] , 
and little research has been performed to expand these models to 
predict genetic merits to account for dominance genetic effects. It 
can be argued that such expansion is difficult because calculation 
becomes complicated and de-regressed estimated breeding values 
are used as phenotypes in most applications of genomic selection 
[9]. However, dominance genetic effect is of theoretical and 
practical important because it is heavily used in crosses of animal 
breeds. In fact, assortative mating and mate allocation boost the 
field performances of livestock [10]. Genomic selection has, 
therefore, renewed the interest in the prediction of dominance 
genetic effects. For example, the dominance genetic variance 
accounted for 5.6% of the phenotypic variance by GBLUP 
including dominance genetic effect [11]. More recendy, GBLUP 
method including dominance genetic effect was suggested and the 



software (GVCBLUP) are already available online (http:// 
animalgene.umn.edu/) [12-14]. 

The present study aimed to evaluate the performance of 
GBLUP including dominance genetic effect by estimating variance 
components and predicting genetic effects for both simulation and 
actual pig data. 

Materials and Methods 

Stochastic Simulation 

A historical population was simulated to establish mutation drift 
equilibrium. The simulated genome comprised one chromosome 1 
Morgan long that contained 6,000 SNP markers and 300 
randomly spaced biallelic quantitative trait loci (QTL). In the 
first generation of the historical population, the initial allele 
frequencies of all markers and QTL were assumed to be 0.5. The 
recurrent mutation process was applied, and the mutation rate of 
markers and QTL was 5.0x10 per locus per generation. 
Recombinations were sampled from a Poisson distribution with a 
mean of 1 per Morgan and were then randomly placed along the 
chromosome. The historical population evolved over 2,000 
generations of random mating and random selection with a 
population size of 100 (50 males and 50 females) to reach 
mutation-drift balance. 

After 2,000 historical generations, a base population (GO) and 
the subsequent 6 generations (Gl to G6) were generated as a 
recent population. The population size of GO increased to 300 
(150 males and 150 females). In Gl to G6, 30 sires were randomly 
selected and mated to 150 dams in each generation. Each dam 
had 1 son and 1 daughter; thus, each sire had 5 sons and 5 
daughters. 
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In GO, 1,000 markers and 50 QTL were randomly selected 
among the segregating markers and QTL with minor allele 
frequencies >0.05. Let Qy and Qy be 2 alleles at the jth QTL. 
The genetic values are then given by a, , 5j, and —a.j for genotypes 
QijQiji QijQy> an d QljQlj, respectively. The value of aj was 
drawn from a gamma distribution with a shape parameter of 0.42; 
its sign was drawn at random with equal chance. According to the 
previous simulation study including dominance genetic effects 
[15], the value of Sj was determined as the product of the absolute 
of «; and the degree of dominance, which was drawn from a 
normal distribution N(0,t 2 ). The total genetic effect (gv) of thejth 
animal was calculated by summing all QTL genotypic values, and 
its variance (ff?) was calculated as the sum of additive and 
dominance genetic variances (a 2 and <Tj) [16], which were 
calculated as follows: 

N q 

<? = ^ 2^(1 -<?,•){«,+ (l-2 qj )dj}\ 



Nq 
7=1 

where qj is the frequency of Qy. The broad-sense heritability (A 2 ) 
of the trait was 0.3. To obtain phenotypic values, an environ- 
mental effect was added to the total genetic effect, which was 
sampled from a normal distribution iV(0,(l — h 2 )a 2 /h 2 ). 

The phenotypes and genotypes of SNP markers were available 
for 1,500 and 1,800 individuals from Gl to G5 and Gl to G6, 
respectively. Thus, the reference population with both phenotypes 
and genotypes comprised 1,500 individuals from Gl to G5, and 
the test population with only genotypes comprised 300 individuals 
in G6. 

In a standard simulation scenario, x and the number of markers 
in GO (Ay were set to 0.5 and 1,000, respectively. To investigate 
the effects of t and N m on the performance of the present method, 
2 alternative scenarios were simulated in addition to the standard 
scenario. Three different values of T (0.25, 0.5, and 1.0) and N m 
(200, 1,000, and 5,000) were simulated in the first and seconds 
groups, respectively. For all of these alternatives, only the intended 
parameter differed from the standard scenario. Twenty replicates 
were simulated for each scenario. 

PIC Pig Data 

Publicly available data including pedigree, genotypic, and 
phenotypic information on a single Pig Improvement Company 
(PIC) nucleus pig line were used (http:/ /www.g3journal.org/ 
content/2/4/429/Suppl/DCl). The total number of individuals 
was 3,512, and all have phenotypes for 2 traits (T4 and T5) and 

Table 1. Summary of PIC pig dataset. 



genotypes available from the PorcineSNP60 chip (jV= 64,233). 
These phenotypes were already adjusted for environmental fixed 
effects such as sex, farm, and year of birth [17]. Then, 1,800 
individuals were randomly selected from all individuals whose 
accuracies from the full PIC dataset exceeded 0.8. From these 
1,800 individuals, 1,500 were selected from old generations and 
defined as the reference population while the other 300 individuals 
were defined as the test population. Genotypes were filtered for 
minor allele frequencies less than 0.05. The pig data are 
summarized in Table 1. 

Genomic BLUP Model 

The GBLUP including additive and dominance genetic effects 
termed GBLUP-D. The statistical model of GBLUP-D can be 

expressed as: 

y = Xb + Za + Zd + e, 

where y is the vector of phenotypes; b is the vector of fixed effects; 
a and d are the vector of additive and dominance genetic effects of 
animals; X and Z are incident matrices for the fixed effects, 
additive, and dominance genetic effect, respectively; and e is the 
vector of residuals. Additive and dominance genetic effects were 
assumed to follow normal distributions: a~A(0,Gc 2 ) and 
d~/V(0,Dff2,), where G and D are additive and dominance 
genomic relationship matrices, respectively. These matrices 
describe the relationships between genotyped individuals and 
can be constructed from the information on genome-wide SNP 
markers. Let Ay and Ay be 2 alleles at thejth marker locus and pj 
be the frequency of Ay. The G matrix is created as follows [3]: 



,,_ M a M a 

J2 2pj(l-Pj) 

7=1 

where M a is the n x N„, matrix (n is a number of individuals) and 
the element of M a for the rth individual at the jth marker is 
calculated as follows: 



r -2 Pj {A x Ax) 
M al/ =J l-2 Pj (AiA 2 ). 
{2-2pj(A 2 A 2 ) 

Similarly, Md is assumed to be the n x /V,„ matrix and the 
element of for the jth individual at the Jth marker can be 
calculated as: 



Accuracy of estimated breeding value 

Number of animals Phenotype from all PIC datasets 

Trait No. SNPs Reference Test Mean SD Mean SD 

T4 27,391 1,500 300 -1.125 2.417 0.875 0.048 

T5 27,287 1,500 300 44.107 60.315 0.880 0.048 

doi:1 0.1 371 /journal.pone.0085792.t001 
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Table 2. Variance component estimates (±standard errors) and heritabilities for simulation data with 3 dominance degrees (0.25, 
0.5, and 1.0). 



Genetic variance components 











Residual 


Narrow-sense 


Broad-sense 


Condition 


Method 


Additive 


Dominance 


variance 


heritability 


heritability 


i = 0.25 


True value 


0.259 


0.032 


0.736 


0.267 


0.299 


N m = 1,000 


GBLUP 


0.277±0.036 




0.793±0.022 


0.256 






GBLUP-D 


0.277±0.036 


0.029±0.012 


0.764±0.024 


0.259 


0.286 


t = 0.5 


True value 


0.191 


0.108 


0.736 


0.185 


0.289 


N m = 1,000 


GBLUP 


0.205±0.029 




0.865±0.024 


0.192 






GBLUP-D 


0.208±0.030 


0.066±0.016 


0.797±0.026 


0.194 


0.256 


T — 1 .0 


True value 


0.138 


0.147 


0.706 


0.139 


0.290 


N m = 1,000 


GBLUP 


0.152±0.024 




0.857±0.024 


0.151 






GBLUP-D 


0.152±0.024 


0.080±0.016 


0.773±0.025 


0.151 


0.231 



doi:1 0.1 371 /journal.pone.0085792.t002 



r -2pj(A x A x ) 

M ihj =\ 2 Pj (\- Pj )(A x A2) . V(d) = M d M d '<T2 = D ( r2, 

[ -2(1 - Pj ) 2 {A 2 A 2 ) 

This element describes the coefficients of 8j in dominance 
deviations. Therefore, d and its variance can be derived as 
follows: 

d = M d 8, 

where 8 is the N m dimensional vector of the jth element, which is 
dj. This dominance formula was also used in GVCBLUP [14]. 
Assuming the dominance genetic effects at different marker loci 
are identically and independently distributed normal variables, the 
variance of the genome-wide dominance effect is calculated as 
follows: 



Table 3. Variance component estimates (±standard errors) and heritabilities for simulation data with 200, 1,000, and 5,000 
markers. 



Genetic variance components 











Residual 


Narrow-sense 


Broad-sense 


Condition 


Method 


Additive 


Dominance 


variance 


heritability 


heritability 


T = 0.5 


True value 


0.195 


0.103 


0.739 


0.188 


0.287 


N m = 200 


GBLUP 


0.158±0.026 




0.878±0.024 


0.153 






GBLUP-D 


0.158±0.026 


0.048±0.012 


0.831 ±0.024 


0.152 


0.199 


t = 0.5 


True value 


0.191 


0.108 


0.736 


0.185 


0.289 


N m = 1,000 


GBLUP 


0.205±0.029 




0.865±0.024 


0.192 






GBLUP-D 


0.208±0.030 


0.066±0.016 


0.797±0.026 


0.194 


0.256 


i = 0.5 


True value 


0.211 


0.075 


0.685 


0.217 


0.295 


N m = 5,000 


GBLUP 


0.206±0.027 




0.761 ±0.021 


0.213 






GBLUP-D 


0.207±0.028 


0.063±0.016 


0.698±0.024 


0.214 


0.279 



doi:1 0.1 371 /journal.pone.0085792.t003 
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7=1 

Consequendy, D can be calculated using M d : 
D= M d M d ' 

N m , ' 

£ {%(!-#)} 

Variance components were estimated with average information 
restricted maximum likelihood (REML) [18]. The dataset of 
reference population were used to predict genetic effects of the 
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Table 4. Accuracies of estimates {r^,, r dd , and r gg ) and 
regression coefficients of estimates on their true values {ba, a , 
b dd , and b gg ) in the test population for simulation data with 3 
dominance degrees (0.25, 0.5, and 1.0). 



Condition 


Method 




r dM 


' g£ 


ba,a 


h,„ 


h« 


i = 0.25 


GBLUP 


0.803 




0.760 


0.898 




0.939 


N m = 1 ,000 


GBLUP-D 


0.804 


0.212 


0.769 


0.902 


0.609 


0.942 


t = 0.5 


GBLUP 


0.743 




0.616 


0.891 




0.976 


N m = 1,000 


GBLUP-D 


0.745 


0.339 


0.664 


0.900 


0.893 


0.994 


T = 1.0 


GBLUP 


0.711 




0.466 


1.001 




0.937 


N m = 1,000 


GBLUP-D 


0.712 


0.478 


0.581 


1.006 


1.189 


1.035 



doi:1 0.1 371 /journal.pone.0085792.t004 



genotyped individuals in test population. The GBLUP-D model 
solutions yielded the estimates of additive (a) and dominance 
genetic effects (d). The estimates of the total genetic effect (g) were 
calculated by the sum of a and d. In GBLUP, g equals a, because 
the dominance genetic effect is not considered. The predictive 
ability of the model was evaluated from accuracy and unbiased- 
ness of estimates in the test population. The accuracies of 
estimated additive, dominance, and total genetic effects (r^ a , 
r dd , and r gg , respectively) were measured as the correlations 
between the estimates and true values. Unbiasedness (ba, a , b dd , 
and bg ig ) was measured using the regressions of estimates on true 
values. A regression coefficient of one denotes unbiasedness. Since 
the true values are unknown in pig data, estimated breeding values 
from the full PIC dataset (a*) and phenotypes iy) were used instead 
of true additive and total genetic effects, respectively. In real data, 
the predictive ability of dominance genetic effect cannot be 
calculated and was inferred from that of total genetic effect. 

Results 

Stochastic Simulation 

Tables 2 and 3 show the estimates of variance components and 
heritability in simulation data with various values of x and N m . For 
all values of X and N„„ there were nearly no differences in 
estimates of additive genetic variance and narrow-sense heritability 
between GBLUP and GBLUP-D. The ratios of dominant genetic 
variance estimated by GBLUP-D were 90.6%, 61.1%, and 54.4% 
with a x of 0.25, 0.5, and 1.0 and 46.6%, 61.1%, and 84.0% with 
an N,„ of 200, 1,000, and 5,000, respectively. 

Tables 4 and 5 show the accuracies and unbiasedness of 
estimated genetic values calculated by GBLUP and GBLUP-D in 
simulation data. For all values of x and N m , r- aa and b^ a were 
almost equal between GBLUP-D and GBLUP. In GBLUP-D, r dd 
and b dd increased with increasing x and N m . Meanwhile, r gtg 
values in GBLUP-D exceeded those in GBLUP by 1.2%, 7.8%, 
and 24.7% with a x of 0.25, 0.5, and 1.0 and by 1.9%, 7.8%, and 
4.2% with an N„, of 200, 1,000, and 5,000, respectively. In 
GBLUP-D, larger values of x resulted in b g g being closer to 1 . 

PIC Pig Data 

In T4 and T5 from the pig data, dominance genetic variance 
accounted for 9.6% and 6.3% of the phenotypic variance, 
respectively (Table 6). The estimated additive genetic variances 
and residual variances calculated by GBLUP-D were smaller than 
those calculated by GBLUP. Thus, GBLUP-D consequently 



Table 5. Accuracies of estimates {r^a, r dd , and r g g ) and 
regression coefficients of estimates on their true values (b~ a a , 
b dd , and b gg ) in the test population for simulation data with 
200, 1,000, and 5,000 markers. 



Condition 


Method 


r a,a 


r d,d 


r g£ 


ba, a 




b is 


t = 0.5 


GBLUP 


0.689 




0.552 


1.049 




1.027 


W m = 200 


GBLUP-D 


0.696 


0.246 


0.563 


1.062 


0.547 


0.938 


t = 0.5 


GBLUP 


0.743 




0.616 


0.891 




0.976 


N m = 1,000 


GBLUP-D 


0.745 


0.339 


0.664 


0.900 


0.893 


0.994 


t = 0.5 


GBLUP 


0.799 




0.694 


1.004 




1.035 


N m = 5,000 


GBLUP-D 


0.801 


0.374 


0.723 


1.010 


0.961 


1.032 



doi:1 0.1 371 /joumal.pone.0085792.t005 



yielded lower narrow-sense heritability and higher broad-sense 
heritability than GBLUP. 

In T4 and T5 from the pig data, there were nearly no 
differences between GBLUP and GBLUP-D with respect to 
r a,a* (Table 7). In T5, r gy and b gy were slightly higher in GBLUP- 
D than GBLUP, whereas r g ^ y and b g<y were not higher in GBLUP- 
D in T4. 

Discussion 

Stochastic Simulation 

The GBLUP-D method captured the substantial ratios of the 
dominance genetic variances and estimated the individual 
dominance genetic effects although there were nearly no 
differences in estimates of additive genetic variance and narrow- 
sense heritability between GBLUP and GBLUP-D. Our result 
indicates that GBLUP-D is expected to improve performance of 
the crossbreds, in particular when degree of dominance is large. 

In the present study, the simulated genome comprised one 
chromosome 1 Morgan long. However, the whole genome sizes of 
livestock are larger. Here, another simulation data were 
constructed to evaluate the effect of the genome size on predictive 
ability. This simulation data comprised five chromosomes of 1 
Morgan. The numbers of markers and QTL set to be 5,000 and 
250 to obtain the same distances between markers and QTL as the 
initial simulation. Table 8 shows the accuracies and unbiasedness 
of estimated genetic values in this simulation data. The dominance 
genetic effects could be captured in this data, but the accuracies of 
additive, dominance and total genetic effects decreased in 
comparison with the genome of 1 chromosome. Hence, large size 
of reference population would be required when genome size is 
large. 

PIC Pig Data 

The degrees of estimated dominance genetic variance in T4 and 
T5 were nearly equal to those in simulation data with a x of 0.25. 
When t was 0.25, r gtg and b g ^ g were 1.2% and 0.3% higher in 
GBLUP-D than GBLUP. These results indicate the predictive 
ability of GBLUP-D in T4 and T5 only improved slightly because 
the degrees of dominance genetic variance were too small. In fact, 
r g y and b g<y in GBLUP-D were little improved in comparison with 
GBLUP. In general, the degree of dominance genetic variance is 
expected to be much larger in crossbred populations than 
purebred ones. Since the present study was based on data from 
purebred PIC pig data, the degrees of dominance genetic variance 
in T4 and T5 might have been small. 
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Table 6. Variance components estimates (± standard errors) and heritabilities for PIC pig data. 





Genetic variance components 


Trait 


Method 


Additive 


Dominance 


Residual variance 


Narrow-sense 
heritability 


Broad-sense 
heritability 


T4 


GBLUP 


1.909 ±0.1 89 




3.678±0.138 


0.342 






GBLUP-D 


1.735±0.201 


0.537±0.219 


3.331 ±0.185 


0.310 


0.405 


T5 


GBLUP 


1298.1 ±123.2 




2198.2±84.3 


0.371 






GBLUP-D 


1239.0± 129.3 


220.8± 125.9 


2049.3±113.1 


0.353 


0.416 


doi:1 0.1 371 /journal.pone.0085792.t006 



Practical Use of GBLUP-D 

GBLUP-D has two practical uses. First, selection on the basis of 
GBLUP-D in crossbreds is useful for commercial production. In 
swine and poultry, crossbreds are the end product. The marker 
information from a purebred and its crossbred relatives enables the 
selection of candidate purebreds for the performance of their 
crossbred offspring [11]. Second, GBLUP-D could allow mating 
allocation to exploit dominance. An extra response will be 
obtained when an appropriate design of future matings using 
mating allocation techniques is implemented [10,19]. 

Additive and Dominance Relationship Matrix 

In the presence of dominance genetic effects, the breeding 
values of Ay Ay, Ay Ay, and Ay Ay at the 7th locus are 
- 2pj{oij + ( 1 - 2 Pj )Sj}, (1 - 2pj){ccj + (1 - 2 Pj )5j} ,and(2 - 2pj) 
{oij + (1 — 2pj)Sj}, respectively. The additive genomic relationship 
matrix in GBLUP-D should be constructed considering <5,-. If all 
QTL are of complete dominance, then dj is 1. However, in 
practice, the value of Sj cannot be determined because the degree 
of dominance is unknown. Therefore, in the present study, the 
additive genomic relationship matrix was the same in GBLUP-D 
and GBLUP. Although there are nearly no differences in the 
estimates of additive genetic effects or their variance between 
GBLUP and GBLUP-D, the additive genomic relationship matrix 
including S may yield good estimates of them. 

A GBLUP method including dominance genetic effects was also 
proposed in a previous study [11]. In that model, the additive 
genomic relationship matrix is same as that in traditional GBLUP. 
However, the dominance genomic relationship matrix (D*) in the 
previous study differs from that in the present study and is defined 
as follows: 



Table 7. Aces of estimates (;"a„»and r^ y ) and regression 
coefficients (ba M > and bg ty ) of a on full PIC dataset {a*) and g 
on phenotypic value (v) in the test population for the PIC pig 
dataset. 



D* 



HH 







'"</,,/■ 




baft' 


b i,y 


Trait 


Method 










T4 


GBLUP 


0.455 


0.286 


0.637 


0.726 




GBLUP-D 


0.456 


0.286 


0.710 


0.724 


T5 


GBLUP 


0.379 


0.288 


0.631 


0.786 




GBLUP-D 


0.376 


0.299 


0.654 


0.805 



£ 2 Pj (l - Pj ){\ -2 Pj (\- Pj )} 
7 = 1 



where the element of H for the zth individual at the jth marker is 
calculated as 



-TpjQ-pjXAtAi) 
11, , = <! \-2 Pj (l -pjXA^). 
-2 Pj {\-p i ){A 2 A 2 ) 



This element corresponds to the heterozygosity coefficient but 
not the dominance genetic effect. To compare the performance of 
GBLUP-D in the present study and the model in the previous 
study [11], predictive ability was calculated in the standard 
simulation scenario (Table 9). GBLUP-D yielded higher accura- 
cies of additive, dominance, and total genetic effects than the 
previous model [11]. This might be because the heterozygosity 
coefficient includes part of the additive genetic effect. In fact, the 
previous study [12] reported that dominance genetic variance 
calculated from the previous model [1 1] was larger than that from 
GBLUP-D. 

Assuming linkage equilibrium and uncorrelated marker effects, 
the dominance genetic variances in the present study {a 2 d ) and the 
previous study [11] (<rjj are calculated as follows: 



Table 8. Accuracies of estimates r dd , and r^ g ) and 
regression coefficients of estimates on their true values (b~ a a , 
b dd , and b gg ) in the test population for simulation data with 3 
dominance degrees (0.25, 0.5, and 1.0) when genome 
comprises of 5 chromosomes with 1 Morgan each. 



doi:1 0.1 371 /journal.pone.0085792.t007 



Condition 


Method 


r aft 


r i4 


r g£ 


ba,a 


h 4 


b is 


I = 0.25 


GBLUP 


0.672 




0.636 


1.012 




1.017 


N m = 1,000 


GBLUP-D 


0.673 


0.148 


0.641 


1.013 


0.810 


1.018 


t = 0.5 


GBLUP 


0.646 




0.502 


1.015 




1.018 


N m = 1,000 


GBLUP-D 


0.647 


0.244 


0.528 


1.016 


1.011 


1.015 


i = 1.0 


GBLUP 


0.612 




0.490 


1.112 




1.124 


N m = 1,000 


GBLUP-D 


0.612 


0.348 


0.525 


1.122 


1.030 


1.120 



doi:1 0.1 371 /journal.pone.0085792.t008 
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Table 9. Accuracies of estimates {r~ a a , r dd , and Yg g ) and 
regression coefficients of estimates on their true values 
(ba, a , b dd , and b^ g ) in the test population for simulation data 
(1 = 0.5,^ = 1,000). 





Method 


I'ti.a 


r J,J 




ba,a 


h,„ 


h i,s 


GBLUP 


0.743 




0.616 


0.891 




0.976 


GBLUP-D 


0.745 


0.339 


0.664 


0.900 


0.893 


0.994 


Su et al. [17] 


0.709 


0.239 


0.655 


0.939 


0.569 


0.967 



doi:1 0.1 371/journal.pone.0085792.t009 
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In addition, the estimated dominance genetic effects in the 
present study (d) and the previous study [11] (d*) are calculated as 
follows: 

d = tT^DZ'V" 1 (y-Xb), 



dW^Z'V-'iy-Xb), 
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where V is the variance of y. If the distribution of the allelic 
frequencies is available, d* can be transformed to d. 

Epistasis 

Increasing knowledge about biological pathways and gene 
networks highlights the importance of gene-gene interactions, i.e., 
epistasis; some authors argue that much of the genetic variance in 
a population is due to such interactions [20-22]. When 
considering second-order epistasis in GBLUP, the epistatic 
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from the Hadamard product of the genomic relationship matrix. 
For example, additive by additive and additive by dominance 
interactions are represented as G#G and G#D, respectively. The 
present study tried to use linear mixed models including G#G and 
G^D for T4 and T5 in pig data. However, the variance 
components of these epistasis could not be estimated because they 
were outside parameter space (data not shown). The previous 
study [11] also included epistatic genomic relationship matrices in 
GBLUP. In that model, the estimated epistatic variances were 
almost 0 for daily gain in Danish Duroc pigs. The estimated 
epistatic variances in the Bayesian model for the percentage of 
CD8 + cells in publicly available mouse data were almost 0 [23]. 
Although these studies could not detect epistastic effect, in recent, 
the marker-generated kinship matrices were suggested in a new 
mixed model method [24] and nonparametric approaches and 
machine-learning techniques were recommended to model more 
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